#!/usr/bin/env python
"""
Draw stars, planets on dbCarta.
Algorithms and stars data are taken from Marble Project at http://edu.kde.org/marble.
Solar System bodies, sattelite's orbit computes by PyEphem module.
"""

__version__ = "0.4"

import calendar, time, math, sys
from __init__ import *
from dbcarta import *
from dbcarta import __doc__ as dd
from demodata.continents import *
from demodata.staricons import *

# PyEphem module for Solar System objects
try:
    import ephem
except:
    print(sys.exc_info()[0], sys.exc_info()[1])

class Qn:
    """Spherical transformations."""
    @staticmethod
    def fromEuler(pitch, yaw, roll):
        cPhi = math.cos(0.5 * pitch)
        cThe = math.cos(0.5 * yaw)
        cPsi = math.cos(0.5 * roll)

        sPhi = math.sin(0.5 * pitch)
        sThe = math.sin(0.5 * yaw)
        sPsi = math.sin(0.5 * roll)

        w = cPhi * cThe * cPsi + sPhi * sThe * sPsi
        x = sPhi * cThe * cPsi - cPhi * sThe * sPsi
        y = cPhi * sThe * cPsi + sPhi * cThe * sPsi
        z = cPhi * cThe * sPsi - sPhi * sThe * cPsi

        return [w, x, y, z]
    @staticmethod
    def toMatrix(qpos):
        vw, vx, vy, vz = qpos
        m = []

        xy = vx * vy
        yy = vy * vy
        zw = vz * vw

        xz = vx * vz
        yw = vy * vw
        zz = vz * vz

        m.append([1.0 - 2.0 * (yy + zz), 2.0 * (xy + zw), 2.0 * (xz - yw), 0.0])

        xx = vx * vx
        xw = vx * vw
        yz = vy * vz

        m.append([2.0 * (xy - zw), 1.0 - 2.0 * (xx + zz), 2.0 * (yz + xw), 0.0])
        m.append([2.0 * (xz + yw), 2.0 * (yz - xw), 1.0 - 2.0 * (xx + yy), 0.0])

        return m
    @staticmethod
    def rotateAroundAxis(centerof, m):
        vw, vx, vy, vz = centerof

        x =  m[0][0] * vx + m[1][0] * vy + m[2][0] * vz
        y =  m[0][1] * vx + m[1][1] * vy + m[2][1] * vz
        z =  m[0][2] * vx + m[1][2] * vy + m[2][2] * vz

        return [1.0, x, y, z]
    @staticmethod
    def inverse(qpos):
        vw, vx, vy, vz = qpos
        inverse = [vw, -vx, -vy, -vz]
        return Qn.normalize(inverse)
    @staticmethod
    def length(qpos):
        vw, vx, vy, vz = qpos
        return math.sqrt(vw * vw + vx * vx + vy * vy + vz * vz)
    @staticmethod
    def normalize(qpos):
        vw, vx, vy, vz = qpos
        v = 1.0 / Qn.length(qpos)
        return [vw * v, vx * v, vy * v, vz * v]
    @staticmethod
    def fromSpherical(lon, lat):
        w = 0.0
        x = math.cos(lat) * math.sin(lon)
        y = math.sin(lat)
        z = math.cos(lat) * math.cos(lon)
        return [w, x, y, z]
    @staticmethod
    def gregorianToJulian(y, m, d):
        if (y <= 99):
            y += 1900
        if (m > 2):
            m -= 3
        else:
            m += 9
            y -= 1
        c = y
        c /= 100
        ya = y - 100*c
        return 1721119 + d + (146097*c)/4 + (1461*ya)/4 + (153*m+2)/5
    @staticmethod
    def siderealTime(gmtm):
        """Greenwich mean sidereal time in hours (0..24)."""
        y, m, d = gmtm[:3]
        mjdUtc = Qn.gregorianToJulian(y, m, d)

        h, m, s = gmtm[3:6]
        offsetUtcSecs = h * 3600 + m * 60 + s
        d_days = mjdUtc - 2451545.5
        d = d_days + ( offsetUtcSecs / ( 24.0 * 3600 ) )

        gmst = 18.697374558 + 24.06570982441908 * d
        return gmst - int( gmst / 24.0 ) * 24.0

class Starry:
    """Render stars, planets, sattelites."""
    starsdata = []
    solar = []
    staricons = []
    gmtm = time.gmtime(time.time()) # UTC date/time set

    def initSky(self, *tags):
        viewport_x, viewport_y = dbcarta.viewportOf()
        self.vx, self.vy = viewport_x * dbcarta.slider.var.get(), viewport_y * dbcarta.slider.var.get()
        self.vx1, self.vy1 = viewport_x * dbcarta.slider['resolution'], viewport_y * dbcarta.slider['resolution']

        self.left, self.right = [self.vx * x for x in dbcarta.dw.xview()]
        self.top, self.bottom = [self.vy * y for y in dbcarta.dw.yview()]

        center_x, center_y = dbcarta.centerOf()[0]
        cx, cy = math.radians(center_x), math.radians(center_y)

        gmst = Qn.siderealTime( self.gmtm )
        skyRotationAngle = gmst / 12.0 * math.pi

        skyAxis = Qn.fromEuler(-cy, cx + skyRotationAngle, 0.0)
        self.skyAxisMatrix = Qn.toMatrix(Qn.inverse(skyAxis))

        self.earthRadius = 180 / math.pi * dbcarta.stepscaleOf() * dbcarta.slider.var.get()
        self.skyRadius = 0.6 * math.sqrt(self.vx1 * self.vx1 + self.vy1 * self.vy1)

        dbcarta.dw.delete(*tags)

    def stars(self, data=[]):
        """Render body on map.
        DATA list of list `[ra, de, mag, nbody, label]`."""
        self.initSky( 'sky' )
        if dbcarta.projectOf() == 2:
            for ra, de, mag, nbody, label in data:
                qpos = Qn.fromSpherical(ra, de)
                w, qx, qy, qz = Qn.rotateAroundAxis(qpos, self.skyAxisMatrix)

                if qz > 0:
                    continue

                earthCenteredX = qx * self.skyRadius
                earthCenteredY = qy * self.skyRadius

                # dark side
                if ( qz < 0 and ( (earthCenteredX * earthCenteredX + earthCenteredY * earthCenteredY) < (self.earthRadius * self.earthRadius) ) ):
                    continue

                x1 = (self.vx1 / 2.0 + self.skyRadius * qx)
                y1 = (self.vy1 / 2.0 - self.skyRadius * qy)

                # outside sky viewport
                if ( (x1 < 0 or x1 >= self.vx1) or (y1 < 0 or y1 >= self.vy1) ):
                    continue

                x = (self.vx / 2.0 + self.skyRadius * qx)
                y = (self.vy / 2.0 - self.skyRadius * qy)

                # outside visible
                if ( (x < self.left or x >= self.right) or (y < self.top or y >= self.bottom) ):
                    continue

                # star point size by magitude
                if ( mag < -1 ): size = 7
                elif ( mag < 0 ): size = 6
                elif ( mag < 1 ): size = 5
                elif ( mag < 2 ): size = 4
                elif ( mag < 3 ): size = 3
                elif ( mag < 4 ): size = 2
                elif ( mag < 5 ): size = 1
                else: size = 0

                if label: text=_(label)
                else: text=None

                dbcarta.dw.create_text([x, y],
                                        anchor='sw',
                                        text=text,
                                        tags=('.body' + str(nbody), 'DotRoute', 'sky'))

                dbcarta.dw.create_image([x, y],
                                        image=self.staricons[size],
                                        tags=('body' + str(nbody), 'DotRoute', 'sky'))

                dbcarta.usercl('body' + str(nbody), {'coords': [['HD',nbody],['ra',ra],['dec',de],['mag',mag]]})

    def sat(self, tledata=[]):
        """Render satellite's orbit (by whorl) on map.
        TLEDATA list of `[sat, line1, line2]`."""
        self.initSky( 'iss' )
        if dbcarta.projectOf() == 2:
            for sat, line1, line2 in tledata:
                iss = ephem.readtle( sat, line1, line2 )
                period = 86400.0 / iss._n # mean motion per sec
                step = 60
                earthRadiusMeters = 6378136
                tcurr = calendar.timegm( self.gmtm )
                tstart = tcurr - period
                tstop = tcurr
                bend = 0
                while not bend:
                    if tstart >= tstop:
                        tstart = tcurr
                        bend = 1

                    t = time.strftime("%Y/%m/%d %H:%M:%S", time.gmtime(tstart))
                    iss.compute(t)
                    xs, ys = math.degrees(float(iss.sublong)), math.degrees(float(iss.sublat))

                    qpos = Qn.fromSpherical(float(iss.ra), float(iss.dec))
                    w, qx, qy, qz = Qn.rotateAroundAxis(qpos, self.skyAxisMatrix)

                    o = ephem.Observer()
                    o.date, o.elevation, o.long, o.lat = t, -earthRadiusMeters, xs, ys
                    iss.compute(o)
                    skyRadius = self.earthRadius * iss.range / earthRadiusMeters

                    tstart += step

                    earthCenteredX = qx * skyRadius
                    earthCenteredY = qy * skyRadius

                    if ( qz < 0 and ( (earthCenteredX * earthCenteredX + earthCenteredY * earthCenteredY) < (self.earthRadius * self.earthRadius) ) ):
                        continue

                    x = (self.vx / 2.0 + skyRadius * qx)
                    y = (self.vy / 2.0 - skyRadius * qy)

                    if ( (x < self.left or x >= self.right) or (y < self.top or y >= self.bottom) ):
                        continue

                    if bend: size = 7
                    else: size = 1

                    if bend:
                        dbcarta.dw.create_image([x, y],
                                                image=self.staricons[size],
                                                tags=('iss' + sat, 'DotRoute', 'iss'))
                        # earth surface
                        dbcarta.loadCarta([('DotPort', 'Surface.' + sat, '((%s,%s))' % (xs, ys), '', '', '', 'yellow')])
                        dbcarta.usercl('iss' + sat, {'coords': [['n',sat],['ra',iss.ra],['dec',iss.dec],['mag',iss.mag]]})
                    else:
                        dbcarta.dw.create_image([x, y],
                                                image=self.staricons[size],
                                                tags=('iss' + sat, 'iss'))

def renderSky():
    """Render sky."""
    # Sun, Moon and planets from PyEphem
    starry.solar = []
    if 'ephem' in sys.modules:
        for body in ( 'Sun', 'Moon', 'Mercury', 'Venus', 'Mars', 'Jupiter', 'Saturn', 'Uranus', 'Neptune', 'Pluto' ):
            o = eval("ephem." + body + "()")
            o.compute('%s/%s/%s %s:%s:%s' % starry.gmtm[:6])
            starry.solar += [[ float(o.ra), float(o.dec), o.mag, body, body ]]
            #print(body, float(o.ra), float(o.dec), o.mag)

        starry.sat( [[ 'ISS (ZARYA)',
                      '1 25544U 98067A   12138.83290882  .00011789  00000-0  17152-3 0  5727',
                      '2 25544  51.6409 274.1861 0010682 328.0828 128.1508 15.56173799773289' ],
                    [ 'MOLNIYA 1-90',
                      '1 24960U 97054A   04104.60038411 -.00000148  00000-0  10000-3 0  8467',
                      '2 24960  64.1699 306.8439 7209074 272.6281  14.3196  2.00635799 48072' ],
                    [ 'COSMOS 158 R/B',
                      '1 02802U 67045B   04105.12224913  .00000082  00000-0  32618-4 0  5785',
                      '2 02802  74.0096 111.2339 0069038 241.3647 118.0552 14.42132256935013' ]] )

    starry.stars( starry.solar + starry.starsdata )
    cx, cy = dbcarta.centerOf()[0]
    centerof = ( math.degrees(P(math.radians(cx), math.pi)),
                 math.degrees(P(math.radians(cy), math.pi/2.0)) * [-1, 1][math.cos(math.radians(cy)) > 0] )
    tmsg_var.set( 'X %s Y %s ' % centerof +
                  'T %s-%s-%s %s:%s:%s' % starry.gmtm[:6] )

def loadStarIcons():
    """Load star size icons."""
    for icon in STARICONS:
        starry.staricons += [BitmapImage(data=icon, foreground='#%02x%02x%02x' % (34, 104, 234))]
    return starry.staricons != []

def loadStars():
    """Load stars data from text."""
    f = open(os.path.join(DEMOPATH, 'demodata', 'stars', 'stars.txt'))
    try:
        for i, buf in enumerate(f):
            ra, de, mag, hip, label = eval(buf)
            starry.starsdata += [[ra, de, mag, hip, label]]
    finally:
        f.close()
    return starry.starsdata != []

def loadConts():
    """Load continents as substrate."""
    data = []
    for label, poly in CONTINENTS:
        for i, p in enumerate(poly):
            type = ['.Mainland', '.Arctic'][label in ('Antarctica', 'Greenland')]
            data += [(type, '%s.%s' % (label,i), p)]
    if data:
        dbcarta.loadCarta(data)
    return data != []

def setTime():
    try:
        msg = [ _('Type date/time in format') + ' "Y-M-D h:m:s":',
              "(" + _('set blank to use current UTC time') + ")" ]
        _gmtm = askstring( _('Set time'), "\n".join(msg), initialvalue=time.strftime("%Y-%m-%d %H:%M:%S", starry.gmtm) )
        if _gmtm:
            _gmtm = time.strptime( _gmtm, "%Y-%m-%d %H:%M:%S" )
        elif _gmtm == None:
            return
        else:
            _gmtm = time.gmtime(time.time())
        starry.gmtm = _gmtm
    except:
        print('setTime: ', sys.exc_info()[0], sys.exc_info()[1])

root = Tk()
dbcarta = dbCarta(root)
_ = setLanguage(dbcarta.langOf(), 'starry')
__ = setLanguage(dbcarta.langOf(), 'dbcarta')
f = Frame(root)
f.pack(fill='both')
Button(f, text=(_("Set time") + "..."), command=lambda : setTime() or renderSky()).pack(side='right')
tmsg_var = StringVar()
Label(f, textvariable=tmsg_var, anchor='w', justify='left').pack(fill='both')
help_var = StringVar()
help_var.set( '\n'.join( [
              _('Turn, Move Starry Sky') + ':',
              '    ' + __('<Shift> + <Arrows> turn sphere (for Globe projection)'),
              '    ' + __('<Shift> + <Left Click> move view with cursor'),
              _('View body info') + ':',
              '    ' + __('<Right Click> context menu')] ) )
Label(root, textvariable=help_var, anchor='w', justify='left').pack(fill='both')
dbcarta.changeProject(2)
dbcarta.loadCarta([('DotPort', 'POV', '((37.61,55.75))', 'POV', '(-59.6423,13.4738)')], docenter=1)
starry = Starry()
if (loadConts() and loadStars() and loadStarIcons()):
    dbcarta.usercl('__master_keyPress', ['renderSky', locals()])
    dbcarta.usercl('__master_keyPressShift', ['renderSky', locals()])
    dbcarta.usercl('__dw_mouseUp', ['renderSky', locals()])
    dbcarta.usercl('__dw_mouseMoveShift', ['renderSky', locals()])
    renderSky()
root.mainloop()
root.destroy()
